Optimizing Urban Logistics with Network Analysis in RotterdamΒΆ
This project explores how network analysis and geospatial data can be used to optimize logistics operations in an urban environment. Using real geospatial datasets of distribution centers in the Netherlands and the road network of Rotterdam, we apply Python libraries such as OSMnx, NetworkX, GeoPandas, and Folium to solve two common logistics problems:
- Shortest Route Optimization between two warehouses
- Closest Facility Assignment for delivery points to their nearest logistics hub
- Service Area Analysis
- Facility Location Optimization
The goal is to simulate real-world supply chain scenarios using data-driven routing strategies and provide interactive visualizations to better understand the results.
The dataset used in this analysis can be found in the following link: https://data.4tu.nl/articles/dataset/Dutch_Distribution_Centres_2021_Geodata/19361018
# Data manipulation
# ==============================================================================
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import plotly.express as px
pd.set_option("display.max_columns", None)
# Data visualisation
# ==============================================================================
import geopandas as gpd
import osmnx as ox
import networkx as nx
import folium
# Load the GeoPackage file and check its structure
# ==============================================================================
gdf = gpd.read_file('distributioncenters2021NL_repo.gpkg')
display(gdf.shape)
display(gdf.columns)
display(gdf.head())
/opt/anaconda3/lib/python3.12/site-packages/pyogrio/geopandas.py:265: UserWarning: More than one layer found in 'distributioncenters2021NL_repo.gpkg': 'buildings2021NL' (default), 'estates2021NL', 'corridor2021NL'. Specify layer parameter to avoid this warning. result = read_func(
(26951, 28)
Index(['osm_id', 'BAG_id', 'Constr_yea', 'Year_cat', 'Footprint', 'Size_cat',
'Function', 'Function_n', 'Job_dens', 'Job_dens_n', 'RIN', 'Year',
'Name', 'Locality', 'Gross_ha', 'Net_ha', 'Max_cat', 'Rail_cat',
'Water_cat', 'Muni_nr', 'Municip', 'Prov_nr', 'Province', 'Corop_nr',
'Corop', 'Corridor', 'Current', 'geometry'],
dtype='object')
| osm_id | BAG_id | Constr_yea | Year_cat | Footprint | Size_cat | Function | Function_n | Job_dens | Job_dens_n | RIN | Year | Name | Locality | Gross_ha | Net_ha | Max_cat | Rail_cat | Water_cat | Muni_nr | Municip | Prov_nr | Province | Corop_nr | Corop | Corridor | Current | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 273364927 | 0441100000002816 | 1972 | 1970s | 778.057432 | S | Trade, Import, Export | 1 | medium | 3 | 1741 | 1971.0 | Witte Paal | Schagen | 76.0000 | 43.0000 | 4 | A | E | 441 | Schagen | 27 | Noord-Holland | 18 | Kop Van Noord-holland | 0 | 0 | MULTIPOLYGON Z (((115780.809 532845.953 0, 115... |
| 1 | 279126121 | 0971100000022139 | 1978 | 1970s | 734.322479 | S | Trade, Import, Export | 1 | medium | 3 | 203678 | 1960.0 | Kerensheide | Stein | 19.8223 | 16.9625 | 3 | D | E | 971 | Stein | 31 | Limburg | 39 | Zuid-limburg | 0 | 1 | MULTIPOLYGON Z (((182791.839 331932.418 0, 182... |
| 2 | 615938117 | 1931100000000715 | 2018 | 2010s and later | 836.482827 | S | Trade, Import, Export | 1 | low | 2 | 171034 | 2012.0 | Zevender | Schoonhoven | 15.0000 | 10.7000 | 4 | A | E | 1931 | Krimpenerwaard | 28 | Zuid-Holland | 28 | Oost-zuid-holland | 0 | 1 | MULTIPOLYGON Z (((118809.492 440320.413 0, 118... |
| 3 | 284074624 | 0302100000119189 | 2013 | 2010s and later | 2048.445693 | S | XXL retail center | 3 | low | 2 | 1203 | 1943.0 | Kom Nunspeet | Nunspeet | 8.7132 | 8.3971 | 5 | A | E | 302 | Nunspeet | 25 | Gelderland | 13 | Veluwe | 0 | 1 | MULTIPOLYGON Z (((181723.773 487021.56 0, 1817... |
| 4 | 280745712 | 0050100000352829 | 1999 | 1990s | 1610.136749 | S | Trade, Import, Export | 1 | very high | 5 | 213 | 1975.0 | Trekkersveld I en II | Zeewolde | 160.1670 | 93.9994 | 5 | A | A | 50 | Zeewolde | 24 | Flevoland | 40 | Flevoland | 0 | 1 | MULTIPOLYGON Z (((163088.239 484983.775 0, 163... |
# Filter the 'gdf' dataset to include 'Transport and Logistic services' and the city of Rotterdam, since it will be the focus of this analysis
# ==============================================================================
rotterdam = gdf[
(gdf['Function'] == 'Transport and logistic services') &
(gdf['Municip'] == 'Rotterdam')
]
# Check new dataframe structure
# ==============================================================================
display(rotterdam.shape)
display(rotterdam.head())
display(rotterdam.tail())
# Define the geographical location to download road network
# ==============================================================================
place_name = 'Rotterdam, Netherlands'
# Download the road network
# ==============================================================================
G = ox.graph_from_place(place_name, network_type= 'drive')
# Plot the street network
# ==============================================================================
ox.plot_graph(G)
(323, 28)
| osm_id | BAG_id | Constr_yea | Year_cat | Footprint | Size_cat | Function | Function_n | Job_dens | Job_dens_n | RIN | Year | Name | Locality | Gross_ha | Net_ha | Max_cat | Rail_cat | Water_cat | Muni_nr | Municip | Prov_nr | Province | Corop_nr | Corop | Corridor | Current | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 33 | 271881098 | 0599100000653535 | 1999 | 1990s | 784.190901 | S | Transport and logistic services | 2 | very high | 5 | 170371 | 2010.0 | Europoort | Rotterdam | 1700.0 | 1369.0 | 5 | B | A | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((72649.646 436919.835 0, 7266... |
| 75 | 273554871 | 0599100000643345 | 1990 | 1990s | 2362.159063 | S | Transport and logistic services | 2 | very high | 5 | 170143 | 1964.0 | Prins Alexander | Rotterdam-noord | 35.0 | 29.0 | 3 | A | E | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((97149.757 440577.129 0, 9718... |
| 181 | 554042637 | None | 1994 | 1990s | 4874.941239 | M | Transport and logistic services | 2 | no data | 99 | 170380 | 1973.0 | Eemhaven | Rotterdam | 358.0 | 292.0 | 4 | B | A | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((88339.665 433065.91 0, 88329... |
| 228 | 273326689 | 0599100010036685 | 1988 | 1980s | 3352.465527 | M | Transport and logistic services | 2 | no data | 99 | 170373 | 1970.0 | Botlek Europoort-oost | Rotterdam | 1200.0 | 1069.0 | 5 | B | A | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((77162.415 434307.927 0, 7718... |
| 229 | 273326686 | 0599100000754059 | 2001 | 2000s | 10457.050191 | L | Transport and logistic services | 2 | low | 2 | 170376 | 1990.0 | Distripark Botlek | Rotterdam | 120.0 | 108.0 | 4 | B | A | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((76763.925 432037.599 0, 7677... |
| osm_id | BAG_id | Constr_yea | Year_cat | Footprint | Size_cat | Function | Function_n | Job_dens | Job_dens_n | RIN | Year | Name | Locality | Gross_ha | Net_ha | Max_cat | Rail_cat | Water_cat | Muni_nr | Municip | Prov_nr | Province | Corop_nr | Corop | Corridor | Current | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 26041 | 273554710 | 0599100000617012 | 1976 | 1970s | 1782.039078 | S | Transport and logistic services | 2 | very low | 1 | 170143 | 1964.0 | Prins Alexander | Rotterdam-noord | 35.0 | 29.0 | 3 | A | E | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((96511.132 440308.944 0, 9654... |
| 26279 | 863262550 | 0599100110003461 | 2020 | 2010s and later | 33920.430577 | XL | Transport and logistic services | 2 | no data | 99 | 170341 | 2008.0 | Maasvlakte | Rotterdam | 2728.0 | 2297.0 | 5 | B | A | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((58428.859 439331.876 0, 5851... |
| 26280 | 863262551 | 0599100110004080 | 2020 | 2010s and later | 20116.635052 | XL | Transport and logistic services | 2 | no data | 99 | 170341 | 2008.0 | Maasvlakte | Rotterdam | 2728.0 | 2297.0 | 5 | B | A | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((58522.069 439072.086 0, 5864... |
| 26445 | 271523822 | 0599100000659860 | 2000 | 2000s | 5113.122114 | M | Transport and logistic services | 2 | very high | 5 | 170164 | 1971.0 | Hoogvliet | Hoogvliet | 28.0 | 19.8 | 3 | A | E | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((85042.696 432049.426 0, 8506... |
| 26896 | 271871679 | 0599100000047089 | 1987 | 1980s | 803.993757 | S | Transport and logistic services | 2 | very high | 5 | 170371 | 2010.0 | Europoort | Rotterdam | 1700.0 | 1369.0 | 5 | B | A | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((71684.971 437914.726 0, 7170... |
(<Figure size 800x800 with 1 Axes>, <Axes: >)
from scipy.spatial import cKDTree
# Convert the graph to GeoDataFrames
# ==============================================================================
nodes, edges = ox.graph_to_gdfs(G)
# Reproject Transport and logistics buildings to match network CRS
# ==============================================================================
rotterdam = rotterdam.to_crs(nodes.crs)
# Build KDTree from network node coordinates
# ==============================================================================
node_coords = np.array([(point.y, point.x) for point in nodes.geometry])
tree = cKDTree(node_coords)
# Get coords of logistics points (using centroids of polygons)
# ==============================================================================
building_coords = np.array([(geom.centroid.y, geom.centroid.x) for geom in rotterdam.geometry])
# Find nearest node index for each index
# ==============================================================================
distances, indices = tree.query(building_coords, k=1)
# Add nearest node ID to the 'rotterdam' dataframe
# ==============================================================================
rotterdam['nearest_node'] = nodes.iloc[indices].index
display(rotterdam.head(3))
display(rotterdam.columns)
| osm_id | BAG_id | Constr_yea | Year_cat | Footprint | Size_cat | Function | Function_n | Job_dens | Job_dens_n | RIN | Year | Name | Locality | Gross_ha | Net_ha | Max_cat | Rail_cat | Water_cat | Muni_nr | Municip | Prov_nr | Province | Corop_nr | Corop | Corridor | Current | geometry | nearest_node | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 33 | 271881098 | 0599100000653535 | 1999 | 1990s | 784.190901 | S | Transport and logistic services | 2 | very high | 5 | 170371 | 2010.0 | Europoort | Rotterdam | 1700.0 | 1369.0 | 5 | B | A | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((4.19026 51.91467 0, 4.19049 ... | 44519439 |
| 75 | 273554871 | 0599100000643345 | 1990 | 1990s | 2362.159063 | S | Transport and logistic services | 2 | very high | 5 | 170143 | 1964.0 | Prins Alexander | Rotterdam-noord | 35.0 | 29.0 | 3 | A | E | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((4.5457 51.95062 0, 4.54624 5... | 44550916 |
| 181 | 554042637 | None | 1994 | 1990s | 4874.941239 | M | Transport and logistic services | 2 | no data | 99 | 170380 | 1973.0 | Eemhaven | Rotterdam | 358.0 | 292.0 | 4 | B | A | 599 | Rotterdam | 28 | Zuid-Holland | 29 | Groot-rijnmond | 1 | 1 | MULTIPOLYGON Z (((4.41902 51.88213 0, 4.41887 ... | 1582380654 |
Index(['osm_id', 'BAG_id', 'Constr_yea', 'Year_cat', 'Footprint', 'Size_cat',
'Function', 'Function_n', 'Job_dens', 'Job_dens_n', 'RIN', 'Year',
'Name', 'Locality', 'Gross_ha', 'Net_ha', 'Max_cat', 'Rail_cat',
'Water_cat', 'Muni_nr', 'Municip', 'Prov_nr', 'Province', 'Corop_nr',
'Corop', 'Corridor', 'Current', 'geometry', 'nearest_node'],
dtype='object')
Shortest Possible RouteΒΆ
In this section, we calculate the shortest path between two randomly selected logistics centers using the Dijkstra algorithm. This is a classic logistics optimization problem focused on minimizing travel distance across a real-world road network.
Key Objectives:
- Learn how to compute optimal routes between fixed locations.
- Understand how to use graph theory (via NetworkX) to traverse road networks.
- Visualize the optimal route on an interactive map using folium.
This analysis is foundational for route planning, helping logistics companies reduce fuel costs and delivery times.
# Define warehouse origin and end point by index (select two random warehouses)
# ==============================================================================
origin_node = rotterdam.iloc[0]['nearest_node']
destination_node = rotterdam.iloc[20]['nearest_node']
# Convert node IDs to lat/lon coordinates
# ==============================================================================
origin_coords = (G.nodes[origin_node]['y'], G.nodes[origin_node]['x'])
destination_coords = (G.nodes[destination_node]['y'], G.nodes[destination_node]['x'])
# Find the shortest path and the shortest path length using NetworkX (apply Dijkstra's algorithm)
# ==============================================================================
shortest_path = nx.shortest_path(G, origin_node, destination_node, weight= 'length')
path_length = nx.shortest_path_length(G, origin_node, destination_node, weight='length')
print(f" Shortest Path: {path_length:.2f} methers")
# Visualise using Folium
# ==============================================================================
m = folium.Map(location=origin_coords, zoom_start=14)
folium.Marker(
location=origin_coords,
icon=folium.Icon(color='blue', icon='info_sign')
).add_to(m) # Marker for origin point
folium.Marker(
location=destination_coords,
icon=folium.Icon(color='green', icon='info_sign')
).add_to(m) # Marker for destination point
route_coords = [(G.nodes[node]['y'], G.nodes[node]['x']) for node in shortest_path]
folium.PolyLine(route_coords, color='red', weight=4, opacity=0.8).add_to(m) # Add the route line
# Display the map
# ==============================================================================
m
Shortest Path: 17558.97 methers
Closest Facility AnalysisΒΆ
In this section, we simulate delivery locations and assign each one to its nearest logistics hub based on road network distance. This is an example of a many-to-one routing problem, commonly used in last-mile delivery planning.
Key Objectives:
- Generate delivery points within the geographic boundary of logistics centers.
- Use nearest neighbor lookup and Dijkstraβs algorithm to determine the optimal warehouse for each point.
- Map the results interactively to visualize how logistics centers serve different delivery zones.
This technique is crucial for logistics companies aiming to balance delivery workloads and reduce response times by choosing the most accessible fulfillment center.
from shapely.geometry import Point
# Get bounding box coordinates of all logistics centers in Rotterdam
# (used to generate random delivery points within the city limits)
# ==============================================================================
minx, miny, maxx, maxy = rotterdam.total_bounds
# Generate 10 random points within the bounding box
# ==============================================================================
num_points = 10
x_coords = np.random.uniform(minx, maxx, num_points)
y_coords = np.random.uniform(miny, maxy, num_points)
delivery_points = [Point(x, y) for x, y in zip(x_coords, y_coords)]
# Create a GeoDataFrame for delivery points
# ==============================================================================
deliveries = gpd.GeoDataFrame(geometry=delivery_points, crs=rotterdam.crs)
# Initialize the map centered around Rotterdam
# ==============================================================================
m = folium.Map(location=[(miny + maxy) / 2, (minx + maxx) / 2], zoom_start=12)
# Build a KDTree of logistics center node coordinates
# ==============================================================================
center_nodes = rotterdam['nearest_node'].values
center_coords = np.array([(G.nodes[node]['x'], G.nodes[node]['y']) for node in center_nodes])
center_tree = cKDTree(center_coords)
# Assign nearest road node to each delivery point
# ==============================================================================
delivery_coords = np.array([(point.x, point.y) for point in deliveries.geometry])
node_ids = list(G.nodes)
node_coords = np.array([(G.nodes[node]['x'], G.nodes[node]['y']) for node in node_ids])
node_tree = cKDTree(node_coords)
# Match each delivery point to its nearest road node
# ==============================================================================
_, node_indices = node_tree.query(delivery_coords, k=1)
deliveries['delivery_node'] = [node_ids[i] for i in node_indices]
# Match each delivery to the closest logistics center based on network distance
# ==============================================================================
delivery_routes = []
for delivery_node in deliveries['delivery_node']:
# Use Dijkstra to all centers
lengths = nx.single_source_dijkstra_path_length(G, delivery_node, weight='length')
# Pick closest center
min_dist = float('inf')
closest_center = None
for center_node in center_nodes:
if center_node in lengths and lengths[center_node] < min_dist:
min_dist = lengths[center_node]
closest_center = center_node
delivery_routes.append({
'delivery_node': delivery_node,
'center_node': closest_center,
'distance': min_dist
})
# Add logistics center to the map
# ==============================================================================
for idx, row in rotterdam.iterrows():
center_coords = (row.geometry.centroid.y, row.geometry.centroid.x)
folium.Marker(location=center_coords, icon=folium.Icon(color='blue', icon='industry')).add_to(m)
# Add delivery points and routes to the map
# ==============================================================================
for idx, row in deliveries.iterrows():
delivery_coords = (row.geometry.y, row.geometry.x)
folium.Marker(location=delivery_coords, icon=folium.Icon(color='green', icon='home')).add_to(m)
# Find the corresponding route
route_info = delivery_routes[idx]
path = nx.shortest_path(G, route_info['delivery_node'], route_info['center_node'], weight='length')
route_coords = [(G.nodes[node]['y'], G.nodes[node]['x']) for node in path]
folium.PolyLine(route_coords, color='red', weight=2.5, opacity=0.8).add_to(m)
# Display the map
# ==============================================================================
m
Service Area AnalysisΒΆ
The Service Area Analysis estimates the geographic coverage or reachability from a logistics center within a certain distance threshold (e.g., 10 km of travel distance). This technique is essential in:
- Territory planning (e.g., determining coverage areas of warehouses)
- Accessibility studies (e.g., identifying underserved areas)
- Network resilience planning (e.g., evaluating which roads are critical to maintain access)
# Choose a few random logistics centers (Origin Nodes)
# ==============================================================================
origin_nodes = rotterdam.iloc[[5, 24, 70]]['nearest_node'].values
colors = ['red', 'green', 'blue']
service_dists = [10000] # or multiple if you want [5000, 10000]
m = folium.Map(location=origin_coords, zoom_start=11)
# Loop over each logistics center to generate and visualize its 10 km service area
# ==================================================================================
for idx, origin_node in enumerate(origin_nodes):
for dist in service_dists:
subgraph = nx.ego_graph(G, origin_node, radius=dist, distance='length')
edges = ox.graph_to_gdfs(subgraph, nodes=False, edges=True)
for _, row in edges.iterrows():
line = folium.PolyLine(
locations=[(y, x) for x, y in row['geometry'].coords],
color=colors[idx],
weight=2.5,
opacity=0.5
)
line.add_to(m)
# Add a marker for the center
# ==============================================================================
origin_coords = (G.nodes[origin_node]['y'], G.nodes[origin_node]['x'])
folium.Marker(location=origin_coords, icon=folium.Icon(color=colors[idx], icon='industry')).add_to(m)
m